Groningen Data - driven haemodynamic response function extraction using Fourier - wavelet regularised

Background: We present a simple, data-driven method to extract haemodynamic response functions (HRF) from functional magnetic resonance imaging (fMRI) time series, based on the Fourier-wavelet regularised deconvolution (ForWaRD) technique. HRF data are required for many fMRI applications, such as defining region-specific HRFs, effciently representing a general HRF, or comparing subject-specific HRFs. Results: ForWaRD is applied to fMRI time signals, after removing low-frequency trends by a wavelet-based method, and the output of ForWaRD is a time series of volumes, containing the HRF in each voxel. Compared to more complex methods, this extraction algorithm requires few assumptions (separability of signal and noise in the frequency and wavelet domains and the general linear model) and it is fast (HRF extraction from a single fMRI data set takes about the same time as spatial resampling). The extraction method is tested on simulated event-related activation signals, contaminated with noise from a time series of real MRI images. An application for HRF data is demonstrated in a simple event-related experiment: data are extracted from a region with significant effects of interest in a first time series. A continuous-time HRF is obtained by fitting a nonlinear function to the discrete HRF coeffcients, and is then used to analyse a later time series. Conclusion: With the parameters used in this paper, the extraction method presented here is very robust to changes in signal properties. Comparison of analyses with fitted HRFs and with a canonical HRF shows that a subject-specific, regional HRF significantly improves detection power. Sensitivity and specificity increase not only in the region from which the HRFs are extracted, but also in other regions of interest. 1 Background In functional magnetic resonance imaging (fMRI), local brain activation temporally changes blood oxygenation, which induces a blood oxygenation level dependent (BOLD) contrast in MR images [1]. Given a model of the BOLD response to a stimulus pattern, statistics can be used to quantify the match between the predicted and measured signals in a voxel, and significant activation is assessed via hypothesis testing. Statistical parametric mapping (SPM) estimates parameters of the noise distribution Published: 10 April 2008 BMC Medical Imaging 2008, 8:7 doi:10.1186/1471-2342-8-7 Received: 17 October 2007 Accepted: 10 April 2008 This article is available from: http://www.biomedcentral.com/1471-2342/8/7 © 2008 Wink et al; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. BMC Medical Imaging 2008, 8:7 http://www.biomedcentral.com/1471-2342/8/7 Page 2 of 16 (page number not for citation purposes) in every voxel to determine a threshold for the computed statistic [2]. It uses the general linear model (GLM): the response to a stimulus pattern is modelled as the output of a linear, time invariant (LTI) system [3]. The response to one type of stimulus is modelled by convolving its temporal distribution with the response function of that type of stimulus. The total response is the sum of responses to all individual stimulus types. The stimulus pattern is known from the experimental setup. However, the haemodynamic response function (HRF), i.e., the temporal BOLD impulse response, is unknown. Essentially, it is a smooth curve that starts to rise two seconds after the stimulus, peaks after six seconds, and returns to baseline within 30 seconds. The resolution of discrete fMRI time samples yields a coarse description of the HRF, which is a problem for designs with short and/or randomised stimuli. In slicewise MRI acquisition, there is no optimal sampling for all slices, and slight differences in temporal onset between voxels in different slices must be accurately modelled to obtain good sensitivity. This requires a continuous-time HRF model. It is sometimes preferable to use a common HRF for many experiments and many regions. In many cases though, a common HRF cannot be assumed [4,5]. Two ways to solve this problem are (i) including a HRF model with more basis functions [5], and (ii) estimating the shape of the HRF in the statistical analysis [4,6]. A drawback of the first and simple approach is that it covers only limited variation in the HRF shape, while more advanced approaches suffer from high complexity and long computation times. For previous work on the variability of the BOLD response between brain areas and even across repeated measurements of the same subject, see Neumann et al. [7]. Extracting a HRF from fMRI data requires assumptions about its shape, and is computationally expensive. A simple method described in the literature is selective averaging with a long interstimulus interval (ISI), assuming nonoverlapping responses [8,9]. Trials with overlapping responses are also averaged, ignoring the fact that the overlaps introduce errors [3,10,11]. Other studies use a functional description of the HRF, whose parameters are determined by curve fitting (examples in [12-14]), one of which [12] uses frequency-domain deconvolution. Another technique based on the GLM expands the HRF into a set of basis functions [15,16]. Ciuciu et al. [4] use a Bayesian method to extract the HRF. The method can simultaneously extract multiple HRFs from multiple experiments. Woolrich et al. [6] use a fully Bayesian approach to fMRI modelling, including the HRF. Both these approaches use certain prior assumptions for the shape of the HRF: in the first case, causality, smoothness, and starting and ending at baseline level, and Gaussian temporal autocorrelations; in the second case, a set of predefined prior HRF shapes, and many priors for modelling the noise distribution and autocorrelation. The HRF extraction method described here is data-driven instead of model-driven. It is based on Fourier-wavelet regularised deconvolution, ForWaRD for short, developed originally for denoising and deblurring of images [17]. ForWaRD combines frequency-domain deconvolution for identifying overlapping signals, frequency-domain regularisation for suppressing noise, and wavelet-domain regularisation for separating signal and noise. It is related to recent wavelet-based deconvolution techniques [18-20], with the advantage that the roles of signal (for fMRI: sparse, high frequency) and response (for fMRI: smooth, low frequency) can be interchanged: unlike other waveletand vaguelette-based deconvolution methods, ForWaRD does the first step entirely in the frequency domain. The method described in this paper uses an fMRI data set and the stimulus time pattern, and produces a time series of image volumes, which contains the HRF in each voxel. Compared to the simple extraction methods above, this method has the advantage of taking overlapping responses into account. On the other hand, it is much simpler than the Bayesian ones described above, in that it does not rely on shape assumptions/priors of the HRF: the extracted time points are determined only by the fMRI signal and the stimuli: the method is data-driven. This important property means that the extracted HRF is not biased by any a priori model. The only requirements for ForWaRD are (i) LTI responses and (ii) separability of signal and noise in the frequency/wavelet domain. Studies of the linearity of the BOLD response indicate that (i) is a reasonable assumption, and given the smoothness of the HRF and the availability of smooth wavelet basis functions, (ii) is satisfied as well. Another important advantage is the fact that wavelet-based signal processing methods are much better in preserving the 1/f-type temporal autocorrelations that are found in fMRI time signals than timedomain methods [21]. Models that do not preserve this structure yield biased results, leading to an increase in false positives [21,22]. Finally, our method is fast: HRF extraction from a single fMRI data set takes about the same time as spatial resampling. To avoid possible misunderstandings we would like to stress that our method only concerns the relation between stimulus and haemodynamic response, but does not make any statement about the underlying neural response. Although it is known from the work by Logothetis et al. [23] that the BOLD response in the visual cortex of the monkey brain roughly corresponds to a convolution of the neuronal local field potential (LFP) BMC Medical Imaging 2008, 8:7 http://www.biomedcentral.com/1471-2342/8/7 Page 3 of 16 (page number not for citation purposes) with some neural response function (NRF), the exact form of this coupling is still unknown and remains an area of active research. In our method, we describe the haemodynamic response as a convolution of the stimulus (not the LFP) with an impulse response function, i.e., the HRF. Given the stimulus pattern our method can recover the HRF, but we do not claim to be able to extract the NRF via deconvolution. The output time series can be used to determine subjectspecific, group-specific, and region-specific HRFs, e.g., by averaging time signals in the corresponding regions. We demonstrate this in an fMRI experiment using two acquisitions, where the HRF extracted from the first fMRI time series is used to predict responses in a second experiment. A functional description for the HRF is found by fitting a combination of two gamma density functions with variable parameters for height, dilation, peak location and lag. We show that contrast and localisation improve by using the extracted HRF instead of the canonical HRF from the SPM2 program (a sum of two gamma density functions (GDFs) with fixed parameters). The remainder of this paper is organised as follows. Section 2.2 reviews the ForWaRD method for regularised deconvolution. Sect. 2.2.4 describes how ForWaRD is used for HRF extraction, and it describes how the method was tested on simulated noisy time series and demonstrates a possible application using event-related fMRI data. Sect. 4 contains some general conclusions. 2 Methods This section briefly reviews the standard General Linear Model (GLM) in fMRI, and describes the ForWaRD method [17] for extracting the HRF. 2.1 The General Linear Model Statistical parametric mapping (SPM) is a common analysis method for fMRI. It (a) computes a statistic for every voxel using the GLM, (b) chooses a threshold based on the parameters of the noise and multiple testing correction; (c) thresholds the statistic map. The GLM describes the response in an fMRI experiment as a weighted sum of explanatory signals. An explanatory signal is the response of one stimulus type. Let the matrix Y[T×N] denote the fMRI data set, whose elements yij have time index i = 1, ..., T and position index j = 1, ..., N . In the GLM,